#!/usr/bin/env perl

=head1 NAME

SubcomplexFinder

=head1 SYNOPSIS

use SubcomplexFinder;
use Scorer;

my $scorer = Scorer->new(maxRMSD=>5, scoreRangeMax=>100);
my $subcomplexFinder = SubcomplexFinder->new(scorer=>$scorer);

my @initialSubcomplexes = $subcomplexFinder(\@queryDomains, \@candidateDomains);

=head1 DESCRIPTION

Identifies subcomplexes of 2 domains based on class identity and similar distances between domains. Receives as input
two lists of domains, each list from a different complex. Domains are represented by the sbglib 'Sphere'-object. Returns 
a list of subcomplexes (subcomplex of two domains).

=head1 REQUIRES

sbglib, PDL

=head1 AUTHOR

Philipp Harth <phil.harth@gmail.com>

=cut

use PDL;
use Scorer;

use Subcomplex;

package SubcomplexFinder;
use Moose;

has 'scorer' => (is => 'rw', isa=> 'Scorer');


=head2 pairwise

 Function: finds subcomplexes of 2 domains  
 Example : 
 Returns : list of subcomplexes 
 Args    : list query domains, list of candidate domains
 
 Tries all non redundant mappings of domain pairs to find subcomplexes of 2 domains based on class identity and similar distance. 
 When the distance deviation of two domain pairs (4 domains, 2 from each complex) is below the maxRMSD threshold, a new subcomplex
 is created. 

=cut
sub pairwise {
	
	#this routine looks for similarity in mutual distances. It would be more efficient, if only domain pairs that are actually in 
	#contact are considered (contact graph instead of distance graph). Contact information could be approximated by using only the 
	#nearest neigbours of a domain.
	
	my @subcomplexes = ();
	
	my ($self, $queryDomains, $candidateDomains)= @_;
	my $nQuery = scalar(@{$queryDomains});
	my $nCandidate = scalar(@{$candidateDomains});
	
	my @queryClasses = map { $_->classification } @$queryDomains;
	my @candidateClasses = map { $_->classification } @$candidateDomains;
	
	my @queryIndexes = ();
	my @candidateIndexes = ();
	
	#test all non-redundant mappings of 2 on 2 domains
	for(my $q1=0; $q1<$nQuery; $q1++){
		for(my $c1=0; $c1<$nCandidate; $c1++){
			next unless $queryClasses[$q1] eq $candidateClasses[$c1]; #additional check required here for wildcard classifications
			#at this point, q1 and c1 are all possible pairs from query and candidate with the same classification
			#start with q2 at $q1+1, because the left side (query) of the mapping can be ordered
			for(my $q2=$q1+1; $q2<$nQuery; $q2++){
				#start with c2 at 0 to cover all symmetric permutations, e.g.: A_1,A_2 <==> A_1',A_2' and A_1,A_2 <==> A_2',A_1'
				for(my $c2=0; $c2<$nCandidate; $c2++){
					next unless $queryClasses[$q2] eq $candidateClasses[$c2]; #additional check required here for wildcard classifications
					if($c1 != $c2){
						
						#test deviation and create subcomplex for mapping
						#in the case of one mapped domain, RMSD is equivalent to euclidean distance
						my $distQuery = SBG::U::RMSD::rmsd($queryDomains->[$q1]->centroid, $queryDomains->[$q2]->centroid);
						my $distCandidate = SBG::U::RMSD::rmsd($candidateDomains->[$c1]->centroid, $candidateDomains->[$c2]->centroid);			
						my $deviation = abs($distQuery-$distCandidate);
		
						if($deviation <= $self->scorer->maxRMSD){
							#create new subcomplex	
							my $score = $self->scorer->scoreRMSD($deviation);
							my $s = Subcomplex->new(numberOfDomains=>2,deviation=>$deviation,rmsd=>$deviation,score=>$score);
							#Set local deviations. Here the domains are not superposed (unlike the domains of the subcomplexes in the merging phase).
							#The two local deviations are therefore set as half of the total deviation.
							my $localDeviations = PDL->new($deviation/2, $deviation/2);
							$s->localDeviations($localDeviations);
							$s->complexA($queryDomains->[0]->pdbid);
							$s->complexB($candidateDomains->[0]->pdbid);
							#add domains to subcomplex, maintaining an alphabetic order based on the domain-ID
							if($queryDomains->[$q1]->id lt $queryDomains->[$q2]->id){
								push @{$s->domainsA()}, $queryDomains->[$q1], $queryDomains->[$q2];
								push @{$s->domainsB()}, $candidateDomains->[$c1], $candidateDomains->[$c2];
							}else{
								push @{$s->domainsA()}, $queryDomains->[$q2], $queryDomains->[$q1];
								push @{$s->domainsB()}, $candidateDomains->[$c2], $candidateDomains->[$c1];
							}
							push @subcomplexes, $s;
						}
						#end of test deviation and create subcomplex
					}
				}
			}
		}
	}
	
	return @subcomplexes;
}

no Moose;
1;

